OHI-Northeast | OHI Science | Citation policy
raw <- read_excel(file.path(dir_anx, "_raw_data/NOAA_NMFS/catch_by_stat_area/Afflerbach_UCSB_Landings by Stat Area w Stock Name & Clam Trips_MAR 2019.xlsx"))
clean <- raw %>%
rename(year = YEAR,
stat_area = `STAT\r\nAREA`,
species = SPECIES,
pounds = `LBS LANDED \r\n(HAIL WT)`,
stock_id = `STOCK ID`,
stock = `STOCK NAME`) %>%
mutate(stat_area = as.numeric(stat_area))
head(clean)## # A tibble: 6 x 6
## year stat_area species pounds stock_id stock
## <dbl> <dbl> <chr> <dbl> <chr> <chr>
## 1 1996 0 CONFIDENTIAL SPECIES 11517161 <NA> <NA>
## 2 1996 462 COD 11827 CODGMSS GOM Cod
## 3 1996 462 CUSK 1065 <NA> <NA>
## 4 1996 462 FLOUNDER, AMERICAN PLAICE… 21191 PLAGMMA Plaice
## 5 1996 462 FLOUNDER, WITCH / GRAY SO… 5496 WITGMMA Witch Floun…
## 6 1996 462 HADDOCK 296 HADGM GOM Haddock
stat_shp <- sf::read_sf(file.path(dir_anx, "spatial/Statistical_Areas_2010_withNames.shp")) %>%
st_set_crs(p4s_nad83) %>%
st_transform(crs = crs(rgns))
stat_shp$stat_area <- st_area(stat_shp) #add area as column
plot(stat_shp["SHORT_NAME"])Overlay stat areas to find what ones are in our NE regions
ohi_stat_areas <- st_intersection(rgns, stat_shp)
ohi_stat_areas$ohi_area <- st_area(ohi_stat_areas)
plot(ohi_stat_areas["Id"])Calculate proportion of each statistical area in our OHI regions
calc_prop_area <- ohi_stat_areas %>%
group_by(Id) %>% #calculate the total statistical area region
mutate(ohi_rgn_prop_area = ohi_area/stat_area) #this column tells us how much of each OHI sub-region falls within the statistical area in our region
plot(calc_prop_area["ohi_rgn_prop_area"])Let’s filter the catch data to just the statistical areas in our region.
region_catch <- clean %>%
filter(stat_area %in% ohi_stat_areas$Id) %>%
left_join(calc_prop_area, by = c("stat_area" = "Id")) %>%
mutate(catch = pounds*ohi_rgn_prop_area) %>%
select(-area_km2, -FULL_NAME, -SHORT_NAME, -stat_area.y, -ohi_area, -NAFODIV)
head(region_catch, n = 20)## # A tibble: 20 x 12
## year stat_area species pounds stock_id stock rgn_name rgn_id state
## <dbl> <dbl> <chr> <dbl> <chr> <chr> <fct> <int> <fct>
## 1 1996 464 COD 31842 CODGMSS GOM … Gulf of… 3 <NA>
## 2 1996 464 CRAB, … 200310 <NA> <NA> Gulf of… 3 <NA>
## 3 1996 464 CUSK 44111 <NA> <NA> Gulf of… 3 <NA>
## 4 1996 464 FLOUND… 13415 PLAGMMA Plai… Gulf of… 3 <NA>
## 5 1996 464 FLOUND… 22237 WITGMMA Witc… Gulf of… 3 <NA>
## 6 1996 464 FLOUND… 2930 YELCCGM CC/G… Gulf of… 3 <NA>
## 7 1996 464 HADDOCK 12570 HADGM GOM … Gulf of… 3 <NA>
## 8 1996 464 HAKE, … 55368 HKSGMNGB Nort… Gulf of… 3 <NA>
## 9 1996 464 HAKE, … 4256 <NA> <NA> Gulf of… 3 <NA>
## 10 1996 464 HAKE, … 25578 HKWGMMA Whit… Gulf of… 3 <NA>
## 11 1996 464 LOBSTE… 439729 <NA> <NA> Gulf of… 3 <NA>
## 12 1996 464 MONKFI… 84788 <NA> <NA> Gulf of… 3 <NA>
## 13 1996 464 POLLOCK 30408 POKGMASS Poll… Gulf of… 3 <NA>
## 14 1996 464 REDFIS… 3960 REDGMGB… Redf… Gulf of… 3 <NA>
## 15 1996 464 SCALLO… 427 <NA> <NA> Gulf of… 3 <NA>
## 16 1996 464 SKATE,… 2104 <NA> <NA> Gulf of… 3 <NA>
## 17 1996 464 WOLFFI… 1481 WOLGMMA Wolf… Gulf of… 3 <NA>
## 18 1996 465 COD 41754 CODGMSS GOM … Gulf of… 3 <NA>
## 19 1996 465 CRAB, … 107500 <NA> <NA> Gulf of… 3 <NA>
## 20 1996 465 CUSK 10134 <NA> <NA> Gulf of… 3 <NA>
## # … with 3 more variables: ohi_rgn_prop_area <dbl>, geometry <POLYGON
## # [m]>, catch <dbl>
Map one species
#black sea bass
bsb <- region_catch %>%
filter(year == 2017, species == "BLACK SEA BASS") %>%
group_by(rgn_id, stock_id, stock, species, year) %>%
summarize(catch = sum(catch))
bsb_map <- rgns_simp %>%
left_join(bsb, by = 'rgn_id')
ggplot(bsb_map) +
geom_sf(aes(fill = catch))+
theme_bw() +
labs(title = "BLACK SEA BASS")Visualize the data
sp_catch <- region_catch %>%
group_by(species, stock, stock_id, year, rgn_name, rgn_id) %>%
summarize(catch = sum(catch))
head(sp_catch, n = 30)## # A tibble: 30 x 7
## # Groups: species, stock, stock_id, year, rgn_name [30]
## species stock stock_id year rgn_name rgn_id catch
## <chr> <chr> <chr> <dbl> <fct> <int> <dbl>
## 1 ALEWIFE <NA> <NA> 1998 Gulf of Maine 3 2627.
## 2 ALEWIFE <NA> <NA> 1998 Maine 6 539.
## 3 ALEWIFE <NA> <NA> 1998 Massachusetts-Gulf of Maine 7 7.07
## 4 ALEWIFE <NA> <NA> 1998 New Hampshire 9 47.8
## 5 ALEWIFE <NA> <NA> 1999 Gulf of Maine 3 43.4
## 6 ALEWIFE <NA> <NA> 1999 Maine 6 8.90
## 7 ALEWIFE <NA> <NA> 1999 Massachusetts-Gulf of Maine 7 0.117
## 8 ALEWIFE <NA> <NA> 1999 New Hampshire 9 0.790
## 9 ALEWIFE <NA> <NA> 2000 Gulf of Maine 3 241.
## 10 ALEWIFE <NA> <NA> 2000 Maine 6 49.4
## # … with 20 more rows
ggplot(sp_catch, aes(x = year, y = catch, color = species)) +
facet_wrap(~rgn_name, scales = "free_y") +
geom_line() +
theme_bw() +
theme(legend.position = 'none')Create maps of catch for all species using just the most recent year
map_catch <- stat_shp %>%
left_join(clean %>%
filter(year == 2017), by = c("Id" = "stat_area")) %>%
filter(!is.na(pounds))
for(i in 1:length(unique(map_catch$species))){
sp <- unique(map_catch$species)[i]
ggplot() +
geom_sf(data = map_catch%>%filter(species == sp), aes(fill = pounds)) +
#geom_sf(data = rgns_simp, aes(), fill = NA) +
theme_bw() +
labs(title = sp)
sp_save_name <- ifelse(str_detect(sp, "/"), str_replace_all(sp, "/", "_"), sp)
ggsave(paste0("figs/", sp_save_name,"_catch_by_statistical_area.pdf"))
}Calculate species catch per OHI region
#first get ohi regions and statistical area proportions
catch_by_ohi_rgn <- region_catch %>%
select(year,species, stock_id, stock, rgn_name, rgn_id, catch) %>%
group_by(species, stock_id, stock, rgn_id, year, rgn_name) %>%
summarize(catch = sum(catch)) %>%
ungroup()
write.csv(catch_by_ohi_rgn, file = "data/nmfs_spatial_catch_by_ohi_rgn.csv")Let’s look at total regional catch for each species (not stock)
#calculate total regional catch per species
species_catch <- catch_by_ohi_rgn %>%
group_by(species, year) %>%
summarize(sp_catch = sum(catch)) %>%
ungroup() %>%
group_by(year) %>%
mutate(yr_catch = sum(sp_catch),
catch_prop = sp_catch/yr_catch) %>%
ungroup() %>%
filter(year > 2004) ggplot(species_catch, aes(x = year, y = catch_prop, fill = species)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(legend.text = element_text(size = 6))There are 155. Way more than we have stock information for. Let’s see what species also have stock ids
stocks <- catch_by_ohi_rgn %>%
filter(!is.na(stock_id)) %>%
select(species, stock, stock_id) %>%
distinct()
stocks## # A tibble: 28 x 3
## species stock stock_id
## <chr> <chr> <chr>
## 1 COD GB Cod East CODGBE
## 2 COD GB Cod West CODGBW
## 3 COD GOM Cod CODGMSS
## 4 CONFIDENTIAL SPECIES Southern Whiting HKSOGMMA
## 5 FLOUNDER, AMERICAN PLAICE /DAB Plaice PLAGMMA
## 6 FLOUNDER, SAND-DAB / WINDOWPANE / BRILL Southern Windowpane FLDSNEMA
## 7 FLOUNDER, SAND-DAB / WINDOWPANE / BRILL Northern Windowpane FLGMGBSS
## 8 FLOUNDER, WINTER / BLACKBACK GB Winter Flounder FLWGB
## 9 FLOUNDER, WINTER / BLACKBACK GOM Winter Flounder FLWGMSS
## 10 FLOUNDER, WINTER / BLACKBACK SNE/MA Winter Flounder FLWSNEMA
## # … with 18 more rows
ggplot(catch_by_ohi_rgn %>%
filter(!is.na(stock_id)), aes(x = year, y = catch, color = stock_id)) +
geom_line() +
facet_wrap(~rgn_name, scales = "free_y") +
theme_bw() +
theme(legend.text = element_text(size = 7))We have records for 0 catch. But what about the missing data? Let’s look at ALEWIFE.
## [1] 1998 1999 2000 2003 2006 2007 2012 2013 2014 2015 2016 2017
Ok clearly we are missing 2001, 2002, 04-05, 2008-11.
We need to gapfill this missing data. When a species/state combination has missing data for a year, we can not assume it has a catch of 0. Since we calculate a rolling average of catch, NAs will remain as NA’s and the average will rely on just one or two years of catch. This is done to account for any wild fluctuations in catch year to year.
gf_data <- catch_by_ohi_rgn %>%
group_by(rgn_id, rgn_name, species, stock, stock_id) %>%
complete(year = 1998:2017) %>%
arrange(year) %>%
mutate(mean_catch = zoo::rollapply(catch, 3, mean, fill = NA, align = 'right')) %>% ## create a new column `mean_catch` with rolling mean of 3 yrs
filter(year > 2004) %>%
select(year, rgn_id, rgn_name, species, stock, stock_id, mean_catch) %>%
ungroup()
ggplot(gf_data, aes(x = year, y = mean_catch, color = species)) +
geom_line() +
facet_wrap(~rgn_name, scales = "free_y") +
theme_bw() +
theme(legend.text = element_text(size = 4),
legend.position = "below")Some of these species are harvested for food as well as other markets like pet food or bait. We want to make sure this goal captures catch meant for human consumption. We have data from NOAA that identifies the amount of catch per species, state and year meant for food, bait, and other markets. This data was cleaned in prop_catch_food_bait.Rmd.
prop_data <- read_csv("data/fish_catch_food_prop_rgn.csv")
toolbox_data <- gf_data %>%
left_join(prop_data) %>%
mutate(prop = ifelse(is.na(prop),1,prop),
mean_catch_times_prop = mean_catch*prop) %>%
filter(!market %in% c("BAIT", "NO MARKET"), #remove bait and no market records.
!is.na(mean_catch)) %>% #remove records with no catch (don't need them)
select(-market, -pounds_live_by_market, -total_pounds_live, -prop)